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FULLY IMPLICIT TIME-STEPPING SCHEMES AND NON-LINEAR SOLVERS FOR 
SYSTEMS OF REACTION-DIFFUSION EQUATIONS 

ANOTIDA MADZVAMUSEt§ AND ANDY H.W. CHUNG* § 


Abstract. In this article we present robust, efficient and accurate fully implicit time-stepping schemes and nonlinear solvers 
for systems of reaction-diffusion equations. The applications of reaction-diffusion systems is abundant in the literature, from 
modelling pattern formation in developmental biology to cancer research, wound healing, tissue and bone regeneration and cell 
motility. Therefore, it is crucial that modellers, analysts and biologists are able to solve accurately and efficiently systems of 
highly nonlinear parabolic partial differential equations on complex stationary and sometimes continuously evolving domains 
and surfaces. The main contribution of our paper is the study of fully implicit schemes by use of the Newton method and the 
Picard iteration applied to the backward Euler, the Crank-Nicolson (and its modifications) and the fractional-step 6 methods. 
Our results conclude that the fractional-step 6 method coupled with a single Newton iteration at each timestep is as accurate 
as the fully adaptive Newton method; and both outperform the Picard iteration. In particular, the results strongly support 
the observation that a single Newton iteration is sufficient to yield as accurate results as those obtained by use of an adaptive 
Newton method. This is particularly advantageous when solving highly complex nonlinear partial differential equations on 
evolving domains and surfaces. To validate our theoretical results, various appropriate numerical experiments are exhibited on 
stationary planary domains and in the bulk of stationary surfaces. 

Key words. Reaction-diffusion systems, fully implicit time-stepping schemes, backward Euler method, Crank-Nicolson 
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1. Introduction Since the pioneering work of Turing JHSj, a wide variety of models of reaction-diffusion 
equations (RDEs) have been proposed as plausible mechanisms for pattern generation processes (191 . Turing 
derived the conditions under which a linearised reaction-diffusion system admits a linearly stable spatially 
homogeneous steady state in the absence of diffusion and yet, in the presence of diffusion, it becomes un¬ 
stable under appropriate conditions to yield a spatially varying inhomogeneous steady state. This process 
is now well-known as diffusively-driven instability and is of particular interest in developmental biological 
pattern formation as a means of initiating self organisation from a virtually homogeneous background. Tur¬ 
ing patterns were first observed by Castets et. al. |7] in a chloride-ionic-malonic-acid reaction. Ouyang 
and Swinney [22 1 were the first to observe a Turing instability from a spatially uniform state to a patterned 
state. Although controversial in a biological context for many years, recent experimental findings suggest this 
may be a mechanism for the formation of repeated structures in skin organ formation mm and zebrafish 
mesoderm cell fates [30] . Beyond developmental biology, reaction-diffusion models are widely applied in cell 
motility, cancer biology, astrophysics, semiconductor physics, ecology, material science, chemistry, financial 
mathematics and textile engineering [TJ 21131 U El Bj 0 [HI EO Q3S EH ■ 

In many cases, these models comprise of highly nonlinear reaction terms which makes it impossible to 
obtain analytical solutions in closed form. Hence numerical methods are employed. For the reaction-diffusion 
theory of pattern formation, these methods play two key roles: (i) they are used to validate linear stability 
theoretical results close to bifurcation points (and vice versa) and (ii) far away from bifurcation points, they 
provide numerical approximate solutions (in the absence of analytical solutions). A typical numerical method 
consists of two processes: first, a space discretisation is applied to render the system of partial differential 
equations (PDEs) into a system of ordinary differentials equations (ODEs) and second, a time discretisation 
is employed, thereby transforming the system of ODEs into a system of linear or nonlinear algebraic equations 
depending on the form of the time discretisation scheme. Finally, techniques from numerical linear algebra 
are employed to solve efficiently the resulting system. 

Space discretisation methods include (but are not limited to) finite differences, finite elements, spectral 
methods, finite volume, closest point methods on stationary domains, volumes and surfaces and more recently 
moving grid, surface finite element and particle methods applied to domains and surfaces that evolve in time 
[mil]. In this article, we choose to use the finite element methodology since it can cope easily with complex 
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irregular geometries, volumes and surfaces. Other numerical methods can be applied, however, for some 
methods (e.g., the finite differences) their extension to complicated, irregular and sometimes continuously 
evolving domains and surfaces is not at all trivial. The applicability of the finite element methods to com¬ 
plicated domains is well known. The finite element method can easily deal with complicated and sometimes 
continuously changing domains. 

Several time discretisation schemes have been widely used to compute solutions of PDEs on stationary 
and evolving domains and surfaces m- These include the forward Euler’s method (the most commonly used 
in computational biology), Gear’s method, a modified Euler predictor-corrector method, Gourlay’s method 
[ 8 ], a semi-implicit Rosenbock integrator [ 8 ], and Runge-Kutta schemes [ 12 ]. Most of these are inadequate 
because of the stiffness of the diffusive term. Fully explicit methods require excessively small time-steps 
resulting in computations that are prohibitively expensive especially in multi-dimensions. Recently, IMEX 
schemes have been used to solve RDEs on stationary one-dimensional domains [TJ [ 23 ] • The key essence 
of these schemes is that an implicit scheme is applied to approximate the diffusive term and an explicit 
scheme is used to approximate the reaction kinetics, hence the name IMEX. More recently, Madzvamuse 
[ 16l presented several time-stepping schemes to compute solutions to reaction-diffusion systems on fixed and 
growing domains. A first order semi-implicit backward Euler differentiation formula (1-SBEM) which treats 
diffusive and linear reaction terms implicitly and nonlinear reaction terms semi-implicitly was introduced 
and shown to be more robust than IMEX schemes. The 1-SBEM employed a single Picard iteration. In all 
these studies, very little work has been done to extend such analysis to fully implicit schemes. Therefore the 
focus of this article is to substantially extend previous IMEX schemes to fully implicit time-stepping schemes 
for RDEs on stationary domains and surfaces. Fully implicit time-stepping schemes offer greater numerical 
stability than the 1-SBEM and IMEX schemes. One revealing result of our studies is that, for RDEs, a single 
Newton iteration is sufficient to obtain as accurate solutions as when an adaptive Newton scheme is used. 
This single Newton iteration for the nonlinear reaction kinetics, outperforms a single Picard iteration applied 
to IMEX time-stepping schemes used in [T 6 ]. On the other hand, the overall elapsed CPU time taken by 
the Picard iteration is substantially lower than that taken by the Newton method for the backward Euler 
and Crank-Nicholson and its modifications. The only exception is the fractional 0-method where the Newton 
method outperforms the Picard iteration in terms of the overall elapsed CPU time. This is attributed to the 
fact that the Newton method uses the GMRES solver to invert the resulting non-symmetric matrices, whereas 
the Picard iteration uses the CG solver to invert block diagonal symmetric matrices. Our recommendation 
is that when solving RDEs on stationary and evolving domains and surfaces using fully implicit schemes, it 
is sufficient to employ single Newton iteration with the fractional 0-method in order to obtain as accurate 
solutions as those obtained when an adaptive Newton method is used. This is particularly relevant for 
problems posed on evolving domains and surfaces ei nn. Numerical experiments demonstrate that the 
fractional 0-method coupled with a single Newton method is about 130 times faster than the backward Euler 
and 15 times faster than the Crank-Nicholson and its modifications. 

Hence, our article is structured as follows: in Section [2] we present the model equations under study 
and the choice of parameter values to be used. The space and time discretisation methods and schemes are 
detailed in Section [3] To validate our theoretical predictions, we present experimental order of convergence 
results in Section [4] Various numerical experiments are presented and discussed in Section [3] In Section 
[ 6 ] we demonstrate how fully implicit schemes and appropriate nonlinear solvers outperform standard IMEX 
schemes. To demonstrate the applicability and generality of the fully implicit scheme, in Section[7]we present 
various solutions of RDEs on stationary planary domains and in the bulk of stationary surfaces. Finally, we 
conclude and outline future research in Section H] 

2. The model equations Let Q be a convex domain with Lipschitz boundary, 1 = [0, T\ be some time 
interval and (x,t) £ D x I. We take for illustrative purposes the well-known Schnakenberg reaction kinetics 
(also known as the activator-depleted substrate or Brusselator model [ 1011251127) 1 to obtain the model system 
of RDEs which reads 


du 

dt 

dv 

dt 


— \ 7 2 u = 7(0 — u + u 2 v) := 7 f(u, v), 

— dS7 2 v = 7 (b — u 2 v) := 7 g(u, v), in fl x I, 


( 2 . 1 ) 


for the pair (u(x,t),v(x,t)) and some real positive numbers a, 5, d and 7 . Here, u and v correspond to 
concentrations of two chemical species, d measures the ratio of the relative diffusivity of the v to u and 7 
measures the strength of the reaction. To close the system, we choose homogeneous Neumann boundary 
conditions on the entire boundary and take initial conditions to be small random perturbations around the 
steady state values 


[Ueq , V e q) — 



(a + b) 2 ) 


( 2 . 2 ) 


Next, we recollect briefly conditions for diffusion-driven instability 1331 used to guide parameter estimation 
as well as yielding linear approximate solutions close to bifurcation points which validate numerical results. 
With no spatial variation in u or v, in the absence of diffusion, the equilibrium point (2.2) is linearly stable 
provided that T9] 


fu 7 fjv 7 0 and fu9v fv9u 7 0 ? 


(2.3) 


where the derivatives are evaluated at the equilibrium point. If one then allows spatial inhomogeneity, it is 
possible that the system evolves to an inhomogeneous steady state. This entails an initial instability caused 
by diffusion, a phenomenon known as diffusion-driven instability, or Turing instability after the author who 
first described it in [33] . Subsequently, the non-linear reaction terms keep the solution bounded in an invariant 
set 2 T)i. 

A linear stability analysis reveals conditions that will drive a Turing instability m ■ Let us consider 
small perturbations from the equilibrium and write them as u{x , t) =u(x : t ) — u eq and v(x, t) = v(x, t) — v eq . 


Define £(x,t) = (u(x,t),v(x,t)). Then linearising (2.1) we obtain the system of linear PDEs 


dt 


= 7 


fu fv 

9u 9v 


1 0 
0 d 


V 2 £, 


(2.4) 


with homogeneous Neumann boundary conditions. This system can be solved by separation of variables to 
yield 


£(*,*) ='^2c k e xt Z k , 


where are the modes which solve the homogeneous Neumann problem 

v 2 a + fc 2 a- = 0. 

These modes will decay with time unless their wavenumber k lies in the range 

k 2 _ < k 2 < fc 2 , 

where 

df u +g v ± \J ( df u + g v ) 2 - 4 d(f u g v - f v g u ) 


k± = 7 ' 

and the following conditions hold: 


2d 


dfu T 9v ^ 0 and ( df u T 4:d(f u g v fvgu') ^ 0? 


(2.5) 

( 2 . 6 ) 

(2.7) 

( 2 . 8 ) 

(2.9) 


where the derivatives are evaluated at the equilibrium values. Thus, if we perturb the system from equilibrium, 
under certain choices of parameters, we can expect exponential growth of some modes which correspond to 
the linearly unstable modes of (2.5). This restriction of parameters defines the so-called Turing space GODDI. 

When Q is the two-dimensional unit square, the eigenmodes of (2.6) have the form cos(n 7 r:r) cos(rmry ) for 
n,m £ Z with eigenvalues k 2 = n 2 + m 2 . The choices a = 0.1, 5 = 0.9, d=10 and 7 = 29 will lead to diffusion- 
driven instability m • The parameters chosen guarantee that the modes corresponding to n 2 + m 2 = 1 are 
linearly unstable. Furthermore, it can be calculated that the corresponding exponential growth rate is about 
1.6246 [19] , This exponential growth rate will later serve as an aid to see if our numerical simulations are 
consistent with the theory. 
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3. Numerical methods Due to the non-linearites, an analytical solution to (2.1 1 is not readily avail¬ 
able, and so we seek numerical solutions using the finite element method as detailed below. 

3.1. Finite element spatial discretisation To discretise in space we use the finite element method 
(FEM) QJG22. We shall experiment with a few time discretisation methods. Time-stepping schemes for this 
reaction system were investigated in m for both stationary and evolving domains. The treatment of the 
non-linear terms was semi-implicit - here we shall treat them implicitly, thereby granting greater stability. 

3.2. Space discretisation Let Th be a triangulation of fi/j C f2, and let each node have coordinates 
Xi, i = 1, 2, • • • , Nh, with Nh denoting the number of degrees of freedom. Let be the set of piecewise 
linear shape functions on Th- Now, implementation of the FEM yields 


{ Mu + jMu + Au — 7 B(u, v)u = yal^, 
Mi) + dAv + 7 B(u , u)v = 7 61^, 

where A and M are the stiffness and mass matrices respectively with entries 


(3.1) 


a ij = / V<^j-V <t>j dx and dx, (3-2) 

J Q J Q 

is the column vector with j-th entry <pj. Given some vectors a and b, B{a 1 b) is the matrix with entries 

(3.3) 


N h N h 

Bij = y ( y &kbi I dx. 

Jn 


fc=1 ;=i 

It is readily checked that given a third vector c, the matrix B satisfies 

B(a, b)c = B(a, c)b = B(c , b)a. 


(3.4) 


Equation (3.1) does not yet lend itself to a numerical solution. First, it is still continuous with respect to 


time. Second, the non-linear matrix B does not allow a solution to be gained in a straightforward manner. 

3.3. Time discretisation For illustrative purposes, we will consider the backward-Euler (BE) and 
Crank-Nicolson (CN) methods and the fractional-step 0-scheme (FSTS) [lTj with a uniform timestep r. To 
proceed, define 


Gi(u, v) = Au + ^ {Mu — 7 B(u, v)u — 7 al^ 
and G 2 (w, v) = dAv + ^B(u, u)v — 76 !^ . 

For the BE discretisation we solve at the (n+l)-th timestep 

? .n+l _ q .n 

M -+ Gi(u n+1 , v n+1 ) = 0, 

T 

v n +1 _ 

M -+ G 2 {u n+ \v n+1 ) = 0. 


(3.5) 

(3.6) 


(3.7) 


For the CN we solve 

„«+1 _ 1 

M ---+ - [Gi(u n+ \v n+1 +Gi(«",»")] 

„.n+l _ ..n | 

M ---+ - [G 2 {u n +\v n+1 ) + G 2 (u n ,v n )\ 

For the FSTS, we split the differential operators, i.e. Gi and G 2 , into two parts. In this problem, there is 
the natural dichotomy of diffusion and reaction; however, here, we shall simply make the distinction between 
terms that are linear and those that are not. Once we have this separation, one splits the timestep into three 


= 0 , 


= 0 . 


(3.8) 
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unequal portions alternately treating each dividend implicitly and explicitly. For some 9 € (0,1/2), we first 
solve for the pair (u n+e , v n+e ) 


M 


u n+0 _ u r, 


0t 


M 


v n+6 _ v n 


Or 


+ Au n+e + -/Mu n+e = 7 al* + 7 B(u n , v n )u n , 
+ dAv n+e = - 7 B(u n , u n )v n , 


(3.9) 


then for the pair ( u n+1 6 , v n+1 9 ) we solve 


M 


M 


u n+l-e _ u n+0 

(1-2 9)t 

v n+1-0 _ v n+0 


-7 B{- 


% n+X-e v „+l-e )u „+l-0 


= 7 al 0 - - 7 Mu n+9 , 


. m + 1 B(u n+1 - e 1 u n+1 - e )v n+1 - e = 7 & 1 * -dAv n+e , 
(1 — 2 6 )t 

and finally for the pair ( u n+1 ,v n+1 ) we solve 


(3.10) 


M 


u n+l _ u n+l- 


M 


_ „.n+l—a 

-h A« n+1 + r )Mu n+1 = yalrf, + 7.B(u n+1_e , ■w rl+1_e )M” +1_6 

dr 

+ dAv" +1 = 7&1 0 - 7 S(M n+1 - e , M n+1 - e )v n+1 - e . 


^n+l _ jjfi+1-9 


(3.11) 


dr 


In the first and third step, the non-linear terms are treated explicitly and the linear terms implicitly. Con¬ 
versely, in the second step, the non-linear terms are treated implicitly and the linear terms explicitly. We 
take 0=1 — 1 /\/2 since it can be shown that this method is second order convergent in time for this value of 

0HH. 

3.4. Techniques for treating the non-linearities To deal with the non-linearities, we shall use the 
Picard iteration and the Newton method [24] . 

At the (n+l)-th timestep, the (Ar-Fl)-th Picard iterate of the BE system (3.71 is the solution of the matrix 
system 


(1+7 )M + A- 7 B« +1 ,< +1 ) 


0 


n+1 
l k +1 


iM + dA + 7BK +1 ,«r 1 ) I l 


'yat rj> + i Mu r ‘ 
7 &l 0 + ^ Mv n 


■ (3.12) 


We note here that the IMEX scheme considered in uni is equivalent to taking a single Picard iteration per 
timestep. Now, take the left-hand side (LHS) of (3.7) and define 


F 1 (u n+1 ,t; n+1 ) = 

F 2 (u n+1 ,v n+1 ) = 


+ 'y)M + A-^B(u n+ \v n+1 ) 


u n+1 — 70 !^- Mu n , 


-M + dA + 7 B(u n+1 ,u n+1 ) 

T 


v n+1 - 7 & 1 U- -Mv n , 

T 


(3.13) 

(3.14) 


and let F = {F 1 ,F 2 ). Then we wish to solve F(u n+1 , v n+1 ) = 0. The (fc + l)-th Newton iterate of the BE 
system is the solution of 




= -F« +1 ,< +1 ), * = 1,2, 


(3.15) 


where Jf is the Jacobian matrix of F and the column vector [u 
be the vertical concatenation of the first and second arguements. Using property (|3.4|) of the matrix B 1 we 


fc+i 1 - u n k +1 ,v n k X\ - < +1 ) is understood to 
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have for some vector £ £ 'R Nh the Gateaux derivative 


Similarly we have, 


&Fi(m 


n+l 

> V k 


du n+1 


-z = 


(l+^M + A-2 1 B(u n + 1 ,v^ 1 ) 


dFAu 


n+l n +1 
k ) v k 


dv n+1 

2 V u k >< +1 ) 


£ = — 7 .B(m 


n+l u n +1 


> “fc 




9F 2 (M+\M n+1 


9m"+ 1 


£ = 2"/B(u 


n+l n+l 


)*> 


and 


9F 2 (m 


n+l n+l 


9n" +1 


* = 


-M + dA + jB(u 


n+l u n+! 


’ “fc 


(3.16) 


(3-17) 

(3.18) 

(3.19) 


Thus, at the (n+l)-th timestep, the (fc+l)-th Newton iterate (3.15) is the solution of the system 

ur\v n k +1 ) - 7 b(«; 


(i+ 7 )M + A-2 7 H( 


n+l _.n+l 




L ) 


U 


2 7 R(m 


n+l -.n+l 


lM + dd + 7 B(«+,+ 1 ) 


n+l 

fc+1 

..n+i 

^fc+1 


/ 7«^ + M M "-2 7 5(+ 1 ,+ 1 K +1 

V 761*+^M«" + 2 7 B« +1 ,< +1 )t;r 1 

We can easily show that the Jacobian is not singular in our case. Indeed, suppose that Jf\i u ™+ 1 v ™+ 1 ) r ) = 0 

for some vector rj£R 2Nh . Writing r/ = (r) 1 , tj 2 ), with rj-,, r 7 2 sK Ar ' 1 , we can add the two resulting equations 
and rearrange to obtain 



Vi = 


1 


+ 7 M + A 


-l 


M + dA I r) 2 


(3.21) 


: = ~K^ 1 K 2 rj 2 , 


where we have used the fact that A and M are both positive definite matrices and are, therefore, invertible. 
Upon substitution, this then implies that 


A7 1 (2 7 13(m+ 1 ,< +1 )AT 1 A' 2 - 7-B« +1 ,< +1 )) 


V2 = V2- 


(3.22) 


If rj 2 += 0, then the matrix on the LHS of (3.22) must be the identity matrix. However, this cannot be 
since the coefficients M+ 1 and v^ +1 appear on the LHS but, obviously, do not appear in the identity matrix. 
Thus, r] 2 must be the zero vector. Similarly, rh must also be the zero vector, from which we conclude that 
Jp is invertible. A similar argument will show that the relevant matrix for Picard iteration is also invertible, 
though it does not indicate whether or not the iterations converge to the desired solution. For this we need 
to show that we have a contraction mapping (for proof see the appendix). 

We can expect linear convergence with the Picard iteration and quadratic convergence with the Newton 
method [22j. The Newton method may fail to converge if the initial guess is far from the solution [Ml . 
However, in our case if we have a small enough timestep, we expect that the solution (u n , v n ) will provide an 
adequate initial guess to find the solution (w" +1 , n" +1 ) since the change in the solution through successive 


timesteps will be small. Evidently, the matrix to be inverted in the Picard system (3.121 is symmetric whereas 


for the Newton system (3.20) it is not. This then precludes the use of efficient iterative techniques such as 
the conjugate gradient method to find a solution to the Newton system. Instead, for the Newton method, 
we use GMRES in order to solve the non-symmetric system of linear algebraic equations [26| . 
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4. Experimental order of convergence (EOC) There is no known analytical solution to the 


Schnakenberg system (2.1 1 . Therefore, we cannot check directly how well our numerical experiments compare 


with the exact solution. In the following, we use a well known method of constructing a solution that will 


satisfy a modified version of ( 2 . 1 ) from which we can indirectly gauge the performance of subsequent discreti¬ 


sation. Though this is no substitute for error estimation, it is, nevertheless, a useful guide. The simulations 
in this section and the next were carried out using the software deal. II [J[. 

In the following, we take the domain Q to be the unit square. Define 




x 

~2 


r 

2 


( 1 + «"*)■ 


(4.1) 


Then u = v = H is the exact solution to the modified system 

^ - V 2 u - 7 (~u + u 2 v) = Hi - V 2 H - y(—H + H 3 ), 


— dS7 2 v + 7 u 2 v = H t — dV 2 ^ -(-7 s 3 , 
at 


(4.2) 


in fl , for t > 0 , 


with homogeneous Neumann boundary conditions. Initial conditions are then defined by 


u 0 = v 0 = 2 ( — - — 


x 3 x 2 


y 3 y 2 


(4.3) 


Here we shall also take the parameter values a = 0.1, & = 0.9, d=10 and 7 = 29 [15]. It is easily seen that u 
and v both tend to the inhomogeneous steady Uq/2 as t—> 00 . 

If we now solve (4.2 1 using the proposed time-step methods we can calculate the error from the exact 
solution H at each timestep. This is done using five different timesteps r, = 2~ l , i = 1,2,3,4, 5 in the time 
interval t £ [0,10]. In the following, the convergence criteria for the Picard iteration and the Newton method 
was chosen to be \\u^\ — m/' + 1 || < 10 -5 , and similarly for the variable v, where ||-|| denotes the L 2 -norm. A 
measure of the error from the exact solution at the (n + l)-th timestep is given by 


i(t 


n+l \ 


n+1 




- E (t 


and a measure of the error of the whole simulation is given by 

E u = e„(10), 


(4.4) 


(4.5) 


which is the error at the end time. Similar quantities can be defined for the ri-variable. We have the error 
estimate [32] 


E u < C[u) (h 2 + r“) , 


(4.6) 


where C[u] is a constant and h is the greatest length of the squares in the mesh T'h- For a particular time¬ 
stepping scheme, the maximal value of a such that the above holds is then its order of convergence. Thus, 
for BE we expect a=l and for the CN and the FSTS we expect a = 2. 

Let E u ,i (* = 1,2,- ,5 ) denote the errors found as above using the same time-stepping scheme with time 

step Tj. The estimate (4.6) is a sharp estimate so that if h 2 ~ r", we have approximately 


OLi 


log(E Uii ) - log(E Uti - 1 ) 
l°g( T i) - log(Ti-l) 


i > 1 , 


(4.7) 


with the aj converging to the expected value as the errors decrease. For the second order methods, the 
coupling hi = Tj was obtained by constructing an n x n square grid with n = 2\ For the backward-Euler 
method hi =, [t\ was taken by constructing a square grid as before but taking the nearest n such that the 
value of 1 /n was closest to yTv 

The errors for r 5 = 1/32 are plotted in Fig. 


4.1 As expected, the second order methods are more 


accurate than BE. The convergence results are shown in Table [4~l| The values of a obtained are in keeping 
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BE 

CN 

FSTS 

OL 1 

-1.965xl0 _s 

2.049 

2.064 

OL2 

1.217 

2.020 

2.020 

OL 3 

0.847 

2.006 

2.006 

CH4 

1.184 

2.002 

2.002 

«5 

0.835 

- 

- 

OLQ 

0.922 

- 

- 

OL'J 

1.083 

- 

- 

ag 

1.048 

- 

- 

OL 9 

0.953 

- 

- 


Table 4.1: Experimental order of convergence (EOC) for the u- variable. Further EOC values were computed 
for BE to investigate the convergence to the expected value of 1. 


with the theory. The apparent non-convergence of the a,; to 1 for BE, due perhaps to the errors not being 
small enough, prompted further values of on to be calculated for BE, keeping the timestep r, = 2 _l . These 
are shown in Table |4~T] where the convergence is made more apparent. 

Upon inspection of Fig. |4.1[ one can observe oscillations in the error for the CN method which are 
damped with time. This behaviour is well known mm , and ultimately mars the results. Put briefly, 
under the CN, high frequency modes of the errors in the initial data do not damp out effectively [2l| . Some 
strategies to reduce these oscillations are studied in and pH] , Here, we use a strategy in which BE is 
performed for a few initial timesteps, afterwhic.h the CN is implemented. Essentially, the BE method damps 
out the oscillations which are otherwise retained in the CN method [6]. We try three variants of this method: 
in the first, we perform BE for the first timestep (CNB); in the second we perform two BE steps (CNB2); 
and in the third we perform five (CNB5). The results of these implementations are shown in Fig. 4.2 for two 
timesteps. These methods seem to have varying results. For r = l/2, CNB seems to increase the oscillations, 
whilst CNB5 reduces them the most, however, for the smaller timestep r = l/32, all methods seems to behave 
similarly. Overall however, the magnitude of the errors are reduced taking more BE timesteps initially. 







(a) (b) 



(c) 


Fig. 4.1: Errors for the modified equation for the u-variable using the adaptive Newton method: (a) BE (b) 
CN (c) FTS. The second order methods CN and FSTS are more accurate than BE. The error ofconvergence 
results are shown in Table 


4.1 
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(a) 


(b) 


Fig. 4.2: Errors in the simulation of the modified equation (4.2) using the modified CN methods CNB, CNB2 
and CNB5 with (a) r=l/2 and (b) r=l/32 (see Section |4| for details). The oscillations observed in CN are 
damped by the modified methods for both these timesteps. An adaptive Newton method was employed for 
the fully implicit time-stepping scheme. 
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5. Numerical Results We have three time-stepping schemes and two ways to solve the non-linearities 
implicitly. We choose Q to be the unit square and use a uniform square mesh with 10,000 squares such that 
the greatest length is h= 1/L00. To ensure errors of a similar magnitude, a timestep of r = 10~ 2 was used 
for the CN method and the FSTS and r = 10 -4 was used for the BE method. For the initial conditions, a 
random perturbation from equilibrium of order ~ 10~ 2 was set on every vertex of the mesh. As before, the 
convergence criterion for the Picard iteration and the Newton method was chosen to be 11— u]! +1 11 < 10~ 5 , 
and similarly for the variable v. 

We expect the system to reach a spatially inhomogeneous steady state. Thus, the simulation was left to 
run until a maximum time f = 30, or until the following criterion were satisfied: 


< 10 


-4 


and 


< 10 


-4 


(5.1) 


This quantity is related to the rate of change of the variables, and so we simply stopped the calculations when 
the solutions have stopped varying significantly with time. We find it convenient to divide by the timestep, 
r, to allow for better comparison between solutions using different timesteps. 

After performing the calculations, we found that the solutions obtained using the Picard iteration and 
the Newton method were virtually indistinguishable. Table |5.1| shows the number of iterations needed to 


reach the steady state as defined by (5.11. For the values quoted for FSTS, the two linear steps which frame 


the non-linear step in each timestep are counted. In Fig. |5.1| the number of non-linear iterations performed 
at each timestep is plotted. The number of Picard iterations varied to a greater degree than those resulting 
from the Newton method. The number of iterations for the Picard iteration varied between 1 and 7, whereas 
for the Newton method it varied between 1 and 2. 




(a) (b) 

Fig. 5.1: Number of non-linear iterations for (a) Picard iteration, n p , and (b) the Newton method, n n , at 
each timestep. 


The convergence history of each time-stepping method using the Newton method is shown in Fig. |5.2| 
The initial growth of the selected mode is evident, as is the decay of the other modes, afterwhich there is 
the evolution to the steady state. This behaviour is reflected in the number of non-linear iterations per 
timestep in Fig. |5.1| as the solution is changing to a greater degree between successive timesteps during the 
exponential growth, more non-linear iterations are needed for convergence. There is good agreement between 
the solutions using the BE method and FSTS, however there is some discrepancy with the results obtained 
using the CN and, in fact, the solution using the CN fails to converge in the chosen time limit. From Fig. 
|5.2| it would seem that there are modes less given to decay in the CN than in the other methods. As seen in 
Section [4j the CN method is prone to oscillations. To address these issues, we shall use the CNB5 method 
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BE 

CN 

CNB5 

FSTS 

end time 

18.45 

30.00 

18.26 

17.91 

total iterations 

201024 

7156 

4353 

3822 

elapsed CPU time (xl0 4 s) 

185.10 

25.92 

10.45 

11.42 


(a) 





BE 

CN 

CNB5 

FSTS 

end time 

18.46 

30.00 

18.48 

18.14 

total iterations 

201052 

5998 

3529 

2994 

elapsed CPU time (xl0 4 s) 

963.6 

105.33 

96.8 

8.39 


(b) 


Table 5.1: Details of the simulation of the Schnakenberg system using (a) Picard iteration and (b) the Newton 
method. Shown are the total number of non-linear iterations required to reach the end time and the elapsed 
CPU time. Note the non-convergence of the CN method. 
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30 


(a) (b) 

Fig. 5.2: Convergence history of the simulation of the Schnakenberg system for (a) the u-variable and (b) the 
v- variable. Seen is the initial decay of the modes, mode excitation (growth) and decay into the inhomogeneous 
steady state. An adaptive Newton method was employed for the fully implicit time-stepping scheme. 


as above where we perform BE for the first five timesteps. The results are shown in Fig. |5.3[ As expected, 
there is much better agreement in this case. 

The theoretical exponential growth rate of the excited mode calculated from the linear stability analysis 
is about A = 1.6246. So, if u ~ e xt we should have 


I u n+1 - u ri 


,n— 11 


0 Ar 


-1 


1 — e 


-At ' 


n > 2. 


(5.2) 


The right-hand side (RHS) is shown in the straight lines in Fig. 5.4 along with the left-hand side (LHS) 


for the variable u in each test. The BE and the FSTS methods show excellent agreement with the theory 
during the period of exponential growth, whereas the CN method shows visible deviation. The results from 
the modified CNB5 method, shown in Fig. 5.^[b) show better agreement. 
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0.2 


FSTS - 



t 

Fig. 5.3: Comparison between the convergence history of the FSTS and the modified CNB5 method for the 
w-variable. Much better agreement is seen here compared to the unmodified CN method. Though not shown, 
similar improvement is seen in the w-variabe as well. An adaptive Newton method was employed for the fully 
implicit time-stepping scheme. 
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(a) 


(b) 


Fig. 5.4: A comparison of growth rates - the horizontal lines correspond to the theoretical growth rate as 
in the RHS of (5.21. (a) BE; (b) CN, CNB5 and FSTS. Excellent agreement during the period of mode 
excitation is seen between the theoretical growth rate and that obtained in the numerical simulations using 
BE, CNB5 and FSTS. Using CN, this agreement is lost. An adaptive Newton method was employed for the 
fully implicit time-stepping scheme. 


6. Comparison with IMEX schemes For the sake of comparison, the simulation was carried out 
without using the iterative techniques to solve the non-linearities but otherwise leaving all other parameters, 
outlined in Section [5j unaltered. Instead, only one Picard iteration was performed at each timestep which is 
equivalent to linearising the non-linear terms to obtain an IMEX scheme. This was done with BE (BIMEX), 
CN (CIMEX), CNB5 (C5IMEX) and FSTS (FIMEX). The results are shown in Table O and Fig. O The 
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Fig. 6.1: 

schemes. 


Comparison of the convergence history of the simulation of the Schnakenberg system for the IMEX 
The solid lines show the IMEX schemes and the dotted lines their fully implicit counterparts. 



BIMEX 

CIMEX 

CBIMEX 

FIMEX 

end time 

18.5024 

30.00 

30.00 

19.15 

elapsed CPU time (xl0 4 s) 

153.7 

16.26 

8.82 

7.65 


Table 6.1: Details of the simulation of the Schnakenberg system using the IMEX schemes with one Picard 
iteration. Shown are the end time of the simulations and the elapsed CPU time. Note the non-convergence 
of CIMEX and C5IMEX. 


results most affected are the CN schemes, both modified and unmodified. Under the IMEX schemes, neither 
converges and both suffer from colossal oscillations spanning magnitudes of order 10 2 . With FIMEX, the 
change is not as dramatic, however in Fig. |6.1| the convergence to the steady state towards the end is visibly 
different. The least affected is BN, with the BIMEX and BN lines in Fig. |6. 1 1 indiscernible from each other. 
This is to be expected since the time-step r is smaller for BN and BIMEX than for the second order methods 
and so the linearisation of the nonlinear terms is more accurate. The solution for CIMEX, C5IMEX and 
FIMEX are compromised due to the fact that being explicit time-stepping schemes, their region of stability is 
reduced [16]. The value r = 0.01 does not fall within this region. However, as we saw in Section [5j this value 
of t gave good results when we treated the non-linearities implicitly. This highlights the fact that in using a 
fully implicit scheme we can allow ourselves to use a larger timestep thereby allowing for greater numerical 
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BE 

CN 

CNB5 

FSTS 

end time 

18.4571 

30.00 

18.49 

18.14 

elapsed CPU time (xl0 4 s) 

793.9 

124.9 

87.3 

9.97 


Table 6.2: Details of the simulation of the Schnakenberg system using the one Newton iteration. Shown are 
the end time of the simulations and the elapsed CPU time. 


stability. 

To further investigate this observation, the simulations were carried out again using only one Newton 
iteration instead of one Picard one. The results agreed well without controlling the number of Newton 
iterations as in Section [5] and the convergence histories are virtually the same. The information about the 
end time of simulation and the CPU time taken is shown in |6.2| Comparing Tables 5.1 and 6.2 


one can 

see good agreement in the end time of simulation. Hence the time taken by the fully adaptive Newton and 
the restricted single Newton methods is similar; there does not seem to be much or any speed up or gained 
by using the fully adaptive versus the single Newton iteration. This is not at all surprising since the full 
Newton varied between one and two iterations per timestep, so only taking one iteration does not actually 
significantly reduce the total elapsed CPU time required. 

7. Applications to other geometries The methods outlined are readily applicable to complex ge¬ 
ometries in multi-dimensions as well as on surfaces. Here we present some solutions of the Schnakenberg 
system in the bulk of the unit sphere, the unit cube and on the unit square for different parameter values. 
For fast and accurate simulations, we use only the FSTS to discretise in time, coupled with the Newton’s 
method (with a single iteration at each timestep) to solve the nonlinearities implicitly. The simulations in 
this section were carried out using the software deal.II [3]. 

We choose parameters to isolate modes on the square and use the same parameters on the 3D geometries. 


As detailed before, for the unit square the eigenmodes of (2.6 1 have the form cos(n7ra) cos^wry) for n,m€Z 


with eigenvalues k 2 = n 2 + m 2 . We choose the parameter values a = 0.1, b = 0.9 and firstly d = 9.1676, 
7=176.72, then secondly d = 8.6076, 7 = 535.09. The first set isolates modes corresponding to (n,m) = (2,1), 
whilst the second isolates those corresponding to (3,3) l5|. The results of such calculations are shown in 
Fig. |7.1| A wide variety of patterns are obtained in these geometries. 
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(e) (0 


Fig. 7.1: Solution for the variable u of the Schnakenberg system on the square (top row), the bulk of the 
cube (middle row) and the bulk of the sphere (bottom row). Parameter values were a = 0.1, 6 = 0.9 and (a), 
(c), (e) d = 9.1676, 7=176.72 and (b), (d), (f) d = 8.6076 ,7 = 535.09. Part of the domain has been cut away 
and shown on the right in the 3D geometries to reveal some internal structure. A single Newton iteration at 
each timestep was employed for the fully implicit time-stepping scheme. 


8. Conclusions Overall elapsed CPU times in numerical experiments shown in Table 


5.1 


demonstrate 


that the fractional 0-method (FSTS) is about 130 times faster than the backward Euler and 15 times faster 
than the Crank-Nicholson and its modifications when both the Picard iteration and the Newton method are 
employed as fully implicit solvers. Furthermore, numerical tests show that a single Newton iteration at each 
timestep outperforms the Picard iteration. However, there is a drawback; from Table [5TT| it is clear that the 
Newton method is slower than the Picard iteration in treating non-linearities arising in reaction-diffusion 
systems when the backward Euler and the Crank-Nicholson method and its modifications are employed. The 
extra time needed is attributed to the need by the Newton method to use the GMRES solver to invert the 
non-symmetric systems, whereas the Picard iteration uses the CG solver for the symmetric systems which is 
quicker. In the FSTS method the non-linear system (3.101 is slightly different than the others and takes less 
iterations for the GMRES solver to solve. 

The BE scheme and FSTS agreed well with each other, however the CN method had to be modified in 
order to gain such an agreement. The modified CN method performs very well; it is second order, is more 
straightforward to implement than the FSTS and is unconditionally stable. However, even though the FSTS 
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does require two extra equations to be solved per timestep, the extra equations are linear and as such do not 
add any significant amount of work to be done to calculate the solution. As a result, we strongly recommend 
using the FSTS coupled with the Newton method (with a single iteration at each timestep) as the preferred 
fully implicit time-stepping scheme for reaction-diffusion systems on stationary and sometimes continuously 
evolving domains, volumes and surfaces. 

The comparison between the IMEX schemes and the fully implicit schemes demonstrate how linearisation 
affects the results. Solving non-linearities implicitly alleviated the need for smaller timesteps. However, this 
alleviation came at the price of solving equations that were non-linear. The choice between IMEX and fully 
implicit schemes then lies in the trade-off between the size of the timestep and the number of extra non-linear 
iterations at each timestep. 

It is to be expected that such results will also apply to other reaction-diffusion systems such as the 
Gierer-Meinhardt model HO], the Thomas model m and Brusselator model H5j- I 11 the case of different 
boundary conditions, the inhomogeneous steady state and the Turing space may change but the same method 
of discretisation may be used, so that the results presented in this article are also applicable. 

Although we have presented results for stationary domains, it is natural to extend the numerical analysis 
of implicit solvers to domains and surfaces that change continuously, i.e. growing domains and evolving 
surfaces. It is when solving non-autonomous systems of PDEs posed on evolving domains and surfaces (as 
well as in the bulk) that a single Newton iteration will become crucially important both as an accurate 
solver as well as an aid for large computational savings. The theory of PDEs posed on evolving domains and 
surfaces (as well as in the bulk) is a fast burgeoning research area with many applications in cell motility, 
plant biology and developmental biology where such models are routinely used 0 OH [31 Ini- 

Appendix. Convergence of the Picard iteration. 

Suppose we have a system of PDE’s of the form 


{ oil , „ 

— - diV u = ci + 7 J[u,v)u, 
dv 

— - d 2 V 2 u = c 2 + 7 <?(«, v)v, in D , for t > 0 , 


(A.l) 


for the diffusion coefficients d\, d 2 > 0, constant source terms Ci, c 2 £ R. and functions / : Q —» R. and 
g: fl—»R. which are in general non-linear, but globally Lipshitz continuous with Lipshitz constants Lf and L g 
respectively, and with no explicit dependence on x or t. Even though the Schnakenberg non-linearities may 
not be globally Lipshitz in general, the solution remains in an invariant set and therefore a Lipschitz constant 
may still be found in this region [341 . We further suppose that u and v always remain positive and that / 
and g are bounded above by some constants /3/ and j3 g respectively. Let Th be a quasi-uniform triangulation 
of H, and Sh be the space of continuous functions on Th that are piecewise linear on each element. Recasting 
into the weak form and using the BE method to discretise time as above we solve for ( u n+1 , v n+1 ) € ShxSh 
from the system 


n+l _ n 

( - ,<t>h) + di(Vu" +1 , Vcfrh) = (ciAh) + (f(u n+ ,v n+1 )u n+1 ,cj) h ), 

T 

n+1 _ ? .71 

( - Ah) + d 2 (Vv n+1 , V(/) h ) = (C2,<t>h) + {g{u n+1 ,V n+1 )v n+l Ah), Vfih € s h , 

T 


(A.2) 


where (•, ■) denotes the inner product in L 2 . We would like to use Picard iteration to solve the non-linear 
equation iteratively, each iteration requiring a linear problem to be solved. The (k + l)-th iterate is then 









+ d l (\/u n k X\A(t>h) 

+ d 2 (V<+\V^) 


(ciAh) + {fAk +1 ’ v k +1 ) u ktlAh), 

Mh) + (gK +1 ,vZ +1 )vEtAh)- 
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(A.3) 








Subtracting (A.2) from (A.3) and choosing appropriate <ph we then have 


~||wfe+l - u n+1 \\' 2 < (f(u k ,Vk)u k + 1 - f(u n+1 ,v n+1 )u n+1 ,u k+1 - u n+1 ) , 
-\\v k +i - u n+1 || 2 < T (g(u k ,v k )v k+ 1 - g(u n+1 ,v n+1 )v n+1 ,u k+1 - u n+1 ) , 


(A.4) 


where || • || denotes the L 2 -norm. We have dropped the superscript (n + 1) from the Picard iterates in (A.4) 
for clarity. Concentrating on the equation in the u-variable, we rewrite the RHS as 


(A.5) 


U k ,v k ) - f(u n+1 ,v n+1 ))(u k+ 1 - u n+1 ),u k+ 1 - u n+1 ^j 
+ 2 [(f(u k ,v k ) - f(u n+1 ,v n+1 ))u n+1 ,u k+ 1 - u n+1 ^ 

+ {f{u n+1 ,v n+1 )u k+ 1 - f(u k , v k )u n+1 , u k+1 - u n+1 ) , 


and obtain 

-|| u fc + 1 - u " +1 || 2 

T 

< ||/(wfc,^fc) - f{u n+1 ,v n+1 )\\(^\\(u k +i - w" +1 ) 2 || + 2||u n+1 || 0O ||u fc +i - u" +1 ||) (A 6) 
+ /3 f \\u k +i - u n+1 || 2 

< L f (Ct + 2|| u "+ 1 || 00 )||u fc +i - u n+1 \\ ||4 - T +1 || + fr|K+i - u n+1 1| 2 , 

for large enough Ci, and where £, k = [u k ,v k ), and similarly for 4~ l ~ 1 - Upon rearrangement, this becomes 

IK+i-«” +1 ||<i/ TZ ^(c 1 + 2|K+ 1 |U)||4-r +1 ll- (a.7) 

Similarly, 

||^ + i-^ +1 ||<A sr ^(C 2 + 2|K +1 || 00 )||^-e ,l+1 ||. (A.8) 

Since ||4+i — 4 +1 || < ||ufc+i — u n+1 || + ||ufe + i — u n+1 || we obtain the estimate 

114+1 - 4 +1 11 < L(c + 2| |4 +1 1 loo) 114 - 4 +1 1|, (A.9) 


where L = rna x{ Lf,L g }, /3 = max{/3y, j3 g }, C = max{Ci, C 2 }. Let 3 be the exact solution (U,V) to the 
original problem ( |A.1| ). Then 

||4 +1 ||oo < ||s(t n+1 ) - 4 +1 ||oo + ||H(i n+1 )||oo < C'(B)(h 2 + t)+ ||E(i" +1 )|| 00> (A.10) 

for some constant C 7 (3). Thus, redefining constants as appropriate, we finally obtain 

114+1 - 4 +1 ll < CL— (1 + C'(T + ft 2 )) 114 - 4 +1 ll 

i-Tp (A. 11 ) 

:=A||4-r +t ||. 

The coefficient A will be less than 1 for small enough r and ft, thus showing that the method converges. In 
particular, A —» 0 as r —> 0 and ft —* 0 simultaneously. 
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